home *** CD-ROM | disk | FTP | other *** search
/ Amiga Plus Leser 15 / Amiga Plus Leser CD 15.iso / Tools / Development / yacas_alg / yacas_morphos / share / yacas / univar.rep / code.ys < prev   
Encoding:
Text File  |  2002-03-13  |  17.4 KB  |  749 lines
  1.  
  2.  
  3.  
  4.  
  5.  
  6. /* todos for univariates:
  7.    - Factorize
  8. */
  9.  
  10.  
  11. /* we need to use "predicates" because we are adding to IsZero */
  12. Use("predicates.rep/code.ys");
  13.  
  14. RuleBase("NormalForm",{expression});
  15. Rule("NormalForm",1,1000,True) expression;
  16.  
  17.  
  18. 0 # NormalForm(UniVariate(_var,_first,_coefs)) <--
  19.     ExpandUniVariate(var,first,coefs);
  20. /*TODO remove
  21. Rule("NormalForm",1,0,Type(expression) = "UniVariate")
  22.     ExpandUniVariate(expression[1],expression[2],expression[3]);
  23. */
  24.  
  25. Function("ExpandUniVariate",{var,first,coefs})
  26. [
  27.   Local(result,i);
  28.   result:=0;
  29.   For(i:=Length(coefs),i>0,i--)
  30.   [
  31.     Local(term);
  32.     term:=NormalForm(coefs[i])*var^(first+i-1);
  33.     result:=result+term;
  34.   ];
  35.   result;
  36. ];
  37.  
  38.  
  39. 10 # IsUniVar(UniVariate(_var,_first,_coefs)) <-- True;
  40. 20 # IsUniVar(_anything) <-- False;
  41.  
  42. /*TODO remove?
  43. Function("IsUniVar",{expr}) Type(expr) = "UniVariate";
  44. */
  45.  
  46. RuleBase("UniVariate",{var,first,coefs});
  47.  
  48. Rule("UniVariate",3,10,Length(coefs)>0 And coefs[1]=0)
  49.   UniVariate(var,first+1,Tail(coefs));
  50. Rule("UniVariate",3,1000,IsComplex(var) Or IsList(var))
  51.     ExpandUniVariate(var,first,coefs);
  52.  
  53. 20 # IsZero(UniVariate(_var,_first,_coefs)) <-- IsZeroVector(coefs);
  54.  
  55. RuleBase("Degree",{expr});
  56. Rule("Degree",1,0, IsUniVar(expr))
  57. [
  58.  
  59.   Local(i,min,max);
  60.   min:=expr[2];
  61.   max:=min+Length(expr[3]);
  62.   i:=max;
  63.   While(i >= min And IsZero(Coef(expr,i))) i--;
  64.   i;
  65. ];
  66.  
  67. 10 # Degree(_poly)      <-- Degree(MakeUni(poly));
  68. 10 # Degree(_poly,_var) <-- Degree(MakeUni(poly,var));
  69.  
  70.  
  71.  
  72. 500 # UniVariate(_var,_f1,_c1) + UniVariate(_var,_f2,_c2) <--
  73. [
  74.   Local(from,result);
  75.   Local(curl,curr,left,right);
  76.  
  77.   Set(curl, f1);
  78.   Set(curr, f2);
  79.   Set(left, c1);
  80.   Set(right, c2);
  81.   Set(result, {});
  82.   Set(from, Min(curl,curr));
  83.  
  84.   While(MathAnd(LessThan(curl,curr),left != {}))
  85.   [
  86.     DestructiveAppend(result,Head(left));
  87.     Set(left,Tail(left));
  88.     Set(curl,MathAdd(curl,1));
  89.   ];
  90.   While(LessThan(curl,curr))
  91.   [
  92.     DestructiveAppend(result,0);
  93.     Set(curl,MathAdd(curl,1));
  94.   ];
  95.   While(MathAnd(LessThan(curr,curl), right != {}))
  96.   [
  97.     DestructiveAppend(result,Head(right));
  98.     Set(right,Tail(right));
  99.     Set(curr,MathAdd(curr,1));
  100.   ];
  101.   While(LessThan(curr,curl))
  102.   [
  103.     DestructiveAppend(result,0);
  104.     Set(curr,MathAdd(curr,1));
  105.   ];
  106.   While(MathAnd(left != {}, right != {}))
  107.   [
  108.     DestructiveAppend(result,Head(left)+Head(right));
  109.     Set(left, Tail(left));
  110.     Set(right, Tail(right));
  111.   ];
  112.   While(left != {})
  113.   [
  114.     DestructiveAppend(result,Head(left));
  115.     Set(left, Tail(left));
  116.   ];
  117.   While(right != {})
  118.   [
  119.     DestructiveAppend(result,Head(right));
  120.     Set(right, Tail(right));
  121.   ];
  122.  
  123.   UniVariate(var,from,result);
  124. ];
  125.  
  126.  
  127. 200 # UniVariate(_var,_first,_coefs) + a_IsNumber <--
  128.       UniVariate(var,first,coefs) + UniVariate(var,0,{a});
  129. 200 # a_IsNumber + UniVariate(_var,_first,_coefs) <--
  130.       UniVariate(var,first,coefs) + UniVariate(var,0,{a});
  131.  
  132. /*TODO remove?
  133. 200 # aLeft_IsUniVar + aRight_IsNumber <-- aRight+aLeft;
  134. 200 # aLeft_IsNumber + aRight_IsUniVar <--
  135.     UniVariate(aRight[1],0,{aLeft})+aRight;
  136. */
  137.  
  138. 200 # - UniVariate(_var,_first,_coefs) <-- UniVariate(var,first,-coefs);
  139.  
  140. /*TODO remove?
  141. 200 # - (aLeft_IsUniVar) <--
  142.      Apply("UniVariate",{aLeft[1],aLeft[2],-(aLeft[3])});
  143. */
  144.  
  145. 200 # aLeft_IsUniVar - aRight_IsUniVar <--
  146. [
  147.   Local(from,result);
  148.   Local(curl,curr,left,right);
  149.  
  150.   curl:=aLeft[2];
  151.   curr:=aRight[2];
  152.   left:=aLeft[3];
  153.   right:=aRight[3];
  154.   result:={};
  155.   from:=Min(curl,curr);
  156.  
  157.   While(curl<curr And left != {})
  158.   [
  159.     DestructiveAppend(result,Head(left));
  160.     left:=Tail(left);
  161.     curl++;
  162.   ];
  163.   While(curl<curr)
  164.   [
  165.     DestructiveAppend(result,0);
  166.     curl++;
  167.   ];
  168.   While(curr<curl And right != {})
  169.   [
  170.     DestructiveAppend(result,-Head(right));
  171.     right:=Tail(right);
  172.     curr++;
  173.   ];
  174.   While(curr<curl)
  175.   [
  176.     DestructiveAppend(result,0);
  177.     curr++;
  178.   ];
  179.   While(left != {} And right != {})
  180.   [
  181.     DestructiveAppend(result,Head(left)-Head(right));
  182.     left  := Tail(left);
  183.     right := Tail(right);
  184.   ];
  185.  
  186.  
  187.   While(left != {})
  188.   [
  189.     DestructiveAppend(result,Head(left));
  190.     left  := Tail(left);
  191.   ];
  192.   While(right != {})
  193.   [
  194.     DestructiveAppend(result,-Head(right));
  195.     right := Tail(right);
  196.   ];
  197.  
  198.   UniVariate(aLeft[1],from,result);
  199. ];
  200.  
  201. /* Repeated squares multiplication
  202.  TODO put somewhere else!!!
  203.  */
  204. 10 # RepeatedSquaresMultiply(_a,- (n_IsInteger)) <-- 1/RepeatedSquaresMultiply(a,n);
  205.  
  206. 15 #  RepeatedSquaresMultiply(UniVariate(_var,_first,{_coef}),(n_IsInteger)) <--
  207.       UniVariate(var,first*n,{coef^n});
  208. 20 # RepeatedSquaresMultiply(_a,n_IsInteger) <--
  209. [
  210.   Local(m,b);
  211.   Set(m,1);
  212.   Set(b,1);
  213.   While(m<=n) Set(m,(ShiftLeft(m,1)));
  214.   Set(m, ShiftRight(m,1));
  215.   While(m>0)
  216.   [
  217.     Set(b,b*b);
  218.     If (MathNot(Equals(BitAnd(m,n), 0)),Set(b,b*a));
  219.     Set(m, ShiftRight(m,1));
  220.   ];
  221.   b;
  222. ];
  223.  
  224. 200 # aLeft_IsUniVar ^ aRight_IsPositiveInteger <--
  225.       RepeatedSquaresMultiply(aLeft,aRight);
  226.  
  227.  
  228.  
  229. /*TODO this can be made twice as fast!*/
  230.  
  231. 200 # (aLeft_IsUniVar * _aRight)_(Not(Contains(VarList(aRight),aLeft[1]))) <--
  232.     aRight*aLeft;
  233.  
  234.  
  235. 200 # (_factor * UniVariate(_var,_first,_coefs))_(Not(Contains(VarList(factor),var))) <--
  236.   UniVariate(var,first,coefs*factor);
  237.  
  238. 200 # (UniVariate(_var,_first,_coefs)/_factor)_(Not(Contains(VarList(factor),var))) <--
  239.   UniVariate(var,first,coefs/factor);
  240.  
  241.  
  242. MaxUniOrder:=3000;
  243. Function("SetOrder",{order}) MaxUniOrder:=order;
  244.  
  245.  
  246. ShiftUniVar(UniVariate(_var,_first,_coefs),_fact,_shift)
  247.    <-- UniVariate(var,first+shift,fact*coefs);
  248.  
  249.  
  250. 200 # UniVariate(_var,_f1,_c1) * UniVariate(_var,_f2,_c2) <--
  251. [
  252.   Local(i,j,n,shifted,result);
  253.   Set(result,MakeUni(0,var));
  254.  
  255.   Set(n,Length(c1));
  256.   For(i:=1,i<=n,i++)
  257.   [
  258.     Set(result,result+ShiftUniVar(UniVariate(var,f2,c2),MathNth(c1,i),f1+i-1));
  259.   ];
  260.   result;
  261. ];
  262.  
  263. /*TODO remove?
  264. Function("ShiftUniVar",{uni,fact,shift})
  265. [
  266.  Apply("UniVariate",{uni[1],uni[2]+shift,fact*(uni[3])});
  267. ];
  268. 200 # (aLeft_IsUniVar * aRight_IsUniVar)_(aLeft[1] = aRight[1]) <--
  269. [
  270.   Local(i,j,n,shifted,result);
  271.   Set(result,MakeUni(0,aLeft[1]));
  272.  
  273.   Set(n,Length(aLeft[3]));
  274.   For(i:=1,i<=n,i++)
  275.   [
  276.     result:=result+ShiftUniVar(aRight,aLeft[3][i],aLeft[2]+i-1);
  277.   ];
  278.   result;
  279. ];
  280. */
  281.  
  282. 5 # Coef(uv_IsUniVar,order_IsList) <--
  283. [
  284.   Local(result);
  285.   result:={};
  286.   ForEach(item,order)
  287.   [
  288.     DestructiveAppend(result,Coef(uv,item));
  289.   ];
  290.   result;
  291. ];
  292.  
  293. 10 # Coef(uv_IsUniVar,_order)_(order<uv[2]) <-- 0;
  294. 10 # Coef(uv_IsUniVar,_order)_(order>=uv[2]+Length(uv[3])) <-- 0;
  295. 20 # Coef(uv_IsUniVar,_order) <-- uv[3][(order-uv[2])+1];
  296. 30 # Coef(uv_CanBeUni,_order) <-- Coef(MakeUni(uv),order);
  297.  
  298. Function("Coef",{expression,var,order})
  299.     NormalForm(Coef(MakeUni(expression,var),order));
  300.  
  301. 10 # LeadingCoef(uv_IsUniVar) <-- Coef(uv,Degree(uv));
  302.  
  303. 20 # LeadingCoef(uv_CanBeUni) <--
  304. [
  305.   Local(uvi);
  306.   uvi:=MakeUni(uv);
  307.   Coef(uvi,Degree(uvi));
  308. ];
  309. 10 # LeadingCoef(uv_CanBeUni,_var) <--
  310. [
  311.   Local(uvi);
  312.   uvi:=MakeUni(uv,var);
  313.   Coef(uvi,var,Degree(uvi));
  314. ];
  315.  
  316.  
  317. Function("UniTaylor",{taylorfunction,taylorvariable,taylorat,taylororder})
  318. [
  319.   Local(n,result,dif,polf);
  320.   result:={};
  321.   [
  322.     MacroLocal(taylorvariable);
  323.     MacroSet(taylorvariable,taylorat);
  324.     DestructiveAppend(result,Eval(taylorfunction));
  325.   ];
  326.   dif:=taylorfunction;
  327.   polf:=(taylorvariable-taylorat);
  328.   For(n:=1,n<=taylororder,n++)
  329.   [
  330.     dif:= Deriv(taylorvariable) dif;
  331.     MacroLocal(taylorvariable);
  332.     MacroSet(taylorvariable,taylorat);
  333.     DestructiveAppend(result,(Eval(dif)/n!));
  334.   ];
  335.   UniVariate(taylorvariable,0,result);
  336. ];
  337.  
  338.  
  339. Function("MakeUni",{expression}) MakeUni(expression,VarList(expression));
  340.  
  341. /* Convert normal form to univariate expression */
  342. RuleBase("MakeUni",{expression,var});
  343.  
  344. 1 # MakeUni(_expr,{}) <-- UniVariate(dummyvar,0,{expression});
  345. 2 # MakeUni(_expr,var_IsList) <-- 
  346. [
  347.   Local(result,item);
  348.   result:=expression;
  349.   ForEach(item,var)
  350.   [
  351.     result:=MakeUni(result,item);
  352.   ];
  353.   result;
  354. ];
  355.  
  356. 10 # MakeUni(UniVariate(_var,_first,_coefs),_var) <--
  357.     UniVariate(var,first,coefs);
  358.  
  359. 20 # MakeUni(UniVariate(_v,_first,_coefs),_var) <--
  360. [
  361.   Local(reslist,item);
  362.   reslist:={};
  363.   ForEach(item,expression[3])
  364.   [
  365.     If(IsFreeOf(item,var),
  366.       DestructiveAppend(reslist,item),
  367.       DestructiveAppend(reslist,MakeUni(item,var))
  368.       );
  369.   ];
  370.   UniVariate(expression[1],expression[2],reslist);
  371. ];
  372.  
  373. /*TODO remove?
  374. Rule("MakeUni",2,10,Type(expression) = "UniVariate")
  375. [
  376.   Local(reslist,item);
  377.   reslist:={};
  378.   ForEach(item,expression[3])
  379.   [
  380.     If(IsFreeOf(item,var),
  381.       DestructiveAppend(reslist,item),
  382.       DestructiveAppend(reslist,MakeUni(item,var))
  383.       );
  384.   ];
  385.   Apply("UniVariate",{expression[1],expression[2],reslist});
  386. ];
  387. */
  388. LocalSymbols(a,b,var,expression)
  389. [
  390.   20 # MakeUni(_expression,_var)_(IsFreeOf(expression,var))
  391.        <-- UniVariate(var,0,{expression});
  392.   30 # MakeUni(_var,_var) <-- UniVariate(var,1,{1});
  393.   30 # MakeUni(_a + _b,_var) <-- MakeUni(a,var) + MakeUni(b,var);
  394.   30 # MakeUni(_a - _b,_var) <-- MakeUni(a,var) - MakeUni(b,var);
  395.   30 # MakeUni(   - _b,_var) <--                - MakeUni(b,var);
  396.   30 # MakeUni(_a * _b,_var) <-- MakeUni(a,var) * MakeUni(b,var);
  397.   30 # MakeUni(_a ^ n_IsInteger,_var) <-- MakeUni(a,var) ^ n;
  398.   30 # MakeUni(_a / _b,_var)_(IsFreeOf(b,var)) <-- MakeUni(a,var) * (1/b);
  399. ];
  400.  
  401. /*TODO remove?
  402. Rule("MakeUni",2,20,IsFreeOf(expression,var)) UniVariate(var,0,{expression});
  403.  
  404. Rule("MakeUni",2,30,expression=var)       UniVariate(var,1,{1});
  405.  
  406. Rule("MakeUni",2,30,Type(expression) = "+")
  407.   MakeUni(expression[1],var)+MakeUni(expression[2],var);
  408. Rule("MakeUni",2,30,Type(expression) = "*")
  409.   MakeUni(expression[1],var)*MakeUni(expression[2],var);
  410. Rule("MakeUni",2,30,Type(expression) = "^" And IsInteger(expression[2]))
  411.   MakeUni(expression[1],var)^expression[2];
  412. Rule("MakeUni",2,30,Type(expression) = "-" And NrArgs(expression) = 1)
  413.      -(MakeUni(expression[1],var));
  414.  
  415. Rule("MakeUni",2,30,Type(expression) = "/" And
  416.      Not(Contains(VarList(expression[2]),var)))
  417.   MakeUni(expression[1],var)*(1/expression[2]);
  418.  
  419. Rule("MakeUni",2,30,Type(expression) = "-" And NrArgs(expression) = 2)
  420.      MakeUni(expression[1],var)-MakeUni(expression[2],var);
  421. */
  422.  
  423.  
  424. 0 # Div(n_IsUniVar,m_IsUniVar)_(Degree(n) < Degree(m)) <-- 0;
  425. 0 # Mod(n_IsUniVar,m_IsUniVar)_(Degree(n) < Degree(m)) <-- n;
  426. 1 # Div(n_IsUniVar,m_IsUniVar)_
  427.     (n[1] = m[1] And Degree(n) >= Degree(m)) <--
  428. [
  429.     UniVariate(n[1],0,
  430.                UniDivide(Concat(ZeroVector(n[2]),n[3]),
  431.                          Concat(ZeroVector(m[2]),m[3]))[1]);
  432. ];
  433. 1 # Mod(n_IsUniVar,m_IsUniVar)_
  434.     (n[1] = m[1] And Degree(n) >= Degree(m)) <--
  435. [
  436.     UniVariate(n[1],0,
  437.                UniDivide(Concat(ZeroVector(n[2]),n[3]),
  438.                          Concat(ZeroVector(m[2]),m[3]))[2]);
  439. ];
  440.  
  441.  
  442.  
  443. /* division algo: (for zero-base univariates:) */
  444. Function("UniDivide",{u,v})
  445. [
  446.   Local(m,n,q,r,k,j);
  447.   m := Length(u)-1;
  448.   n := Length(v)-1;
  449.   While (m>0 And IsZero(u[m+1])) m--;
  450.   While (n>0 And IsZero(v[n+1])) n--;
  451.   q := ZeroVector(m-n+1);
  452.   r := FlatCopy(u);  /*  (m should be >= n) */
  453.   For(k:=m-n,k>=0,k--)
  454.   [
  455.     q[k+1] := r[n+k+1]/v[n+1];
  456.     For (j:=n+k-1,j>=k,j--)
  457.     [
  458.       r[j+1] := r[j+1] - q[k+1]*v[j-k+1];
  459.     ];
  460.   ];
  461.   Local(end);
  462.   end:=Length(r);
  463.   While (end>n)
  464.   [
  465.     DestructiveDelete(r,end);
  466.     end:=end-1;
  467.   ];
  468.  
  469.   {q,r};
  470. ];
  471.  
  472.  
  473. DropEndZeroes(list):=
  474. [
  475.   Local(end);
  476.   end:=Length(list);
  477.   While(list[end] = 0)
  478.   [
  479.     DestructiveDelete(list,end);
  480.     end:=end-1;
  481.   ];
  482. ];
  483.  
  484.  
  485.  
  486. Function("UniGcd",{u,v})
  487. [
  488.   Local(l,div,mod,m);
  489.  
  490.   DropEndZeroes(u);
  491.   DropEndZeroes(v);
  492. /*
  493.   If(Length(v)>Length(u),
  494.     [
  495.       Locap(swap);
  496.       swap:=u;
  497.       u:=v;
  498.       v:=swap;
  499.     ] );
  500.   If(Length(u)=Length(v) And v[Length(v)] > u[Length(u)],
  501.     [
  502.       Locap(swap);
  503.       swap:=u;
  504.       u:=v;
  505.       v:=swap;
  506.     ] );
  507.   */
  508.  
  509.  
  510.   l:=UniDivide(u,v);
  511.  
  512.   div:=l[1];
  513.   mod:=l[2];
  514.  
  515.   DropEndZeroes(mod);
  516.   m := Length(mod);
  517.  
  518. /* Echo({"v,mod = ",v,mod}); */
  519. /*  If(m <= 1, */
  520.   If(m = 0,
  521.      v,
  522. /*     v/v[Length(v)], */
  523.      UniGcd(v,mod));
  524. ];
  525.  
  526.  
  527.  
  528. 0 # Gcd(n_IsUniVar,m_IsUniVar)_
  529.     (n[1] = m[1] And Degree(n) < Degree(m)) <-- Gcd(m,n);
  530.  
  531. 1 # Gcd(nn_IsUniVar,mm_IsUniVar)_
  532.     (nn[1] = mm[1] And Degree(nn) >= Degree(mm)) <--
  533. [
  534.    UniVariate(nn[1],0,
  535.                 UniGcd(Concat(ZeroVector(nn[2]),nn[3]),
  536.                        Concat(ZeroVector(mm[2]),mm[3])));
  537. ];
  538.  
  539. RuleBase("PSolve",{uni});
  540.  
  541. Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 1)
  542.     -Coef(uni,0)/Coef(uni,1);
  543.  
  544. Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 2)
  545.     [
  546.      Local(a,b,c,d);
  547.      c:=Coef(uni,0);
  548.      b:=Coef(uni,1);
  549.      a:=Coef(uni,2);
  550.      d:=b*b-4*a*c;
  551.      {(-b+Sqrt(d))/(2*a),(-b-Sqrt(d))/(2*a)};
  552.     ];
  553.  
  554.  
  555. Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 3 )
  556.     [
  557.      Local(p,q,r,w,ww,a,b);
  558.      Local(coef0,coef1,coef3,adjust);
  559.  
  560. /* Get coefficients for a new polynomial, such that the coefficient of
  561.    degree 2 is zero:
  562.    Take f(x)=a0+a1*x+a2*x^2+a3*x^3 and substitute x = x' + adjust
  563.    This gives g(x) = b0+b1*x+b2*x^2+b3*x^3 where
  564.    b3 = a3;
  565.    b2 = 0 => adjust = (-a2)/(3*a3);
  566.    b1 = 2*a2*adjust+3*a3*adjust^2+a1;
  567.    b0 = a2*adjust^2+a3*adjust^3+adjust*a1+a0;
  568.  
  569.    After solving g(x') = 0, return x = x' + adjust.
  570. */
  571.  
  572.      adjust := (-Coef(uni,2))/(3*Coef(uni,3));
  573.      coef3 := Coef(uni,3);
  574.      coef1 := 2*Coef(uni,2)*adjust+3*Coef(uni,3)*adjust^2+Coef(uni,1);
  575.      coef0 := Coef(uni,2)*adjust^2+Coef(uni,3)*adjust^3+
  576.               adjust*Coef(uni,1)+Coef(uni,0);
  577.  
  578.      p:=coef3;
  579.      q:=coef1/p;
  580.      r:=coef0/p;
  581.     w:=Complex(-1/2,Sqrt(3/4));
  582.     ww:=Complex(-1/2,-Sqrt(3/4));
  583.  
  584. /* Equation is xxx + qx + r = 0 */
  585. /*     Let x = a + b
  586.     a^3 + b^3 + 3(aab + bba) + q(a + b) + r = 0
  587.     a^3 + b^3 + (3ab+q)x + r = 0
  588.  
  589.     Let 3ab+q = 0. This is permissible, for we can still find a+b == x
  590.  
  591.     a^3 + b^3 = -r
  592.     (ab)^3 = -q^3/27
  593.  
  594.         So a^3 and b^3 are the roots of t^2 + rt - q^3/27 = 0
  595.  
  596.     Let
  597.                 a^3 = -r/2 + Sqrt(q^3/27+ rr/4)
  598.                 b^3 = -r/2 - Sqrt(q^3/27+ rr/4)
  599.     Therefore there are three values for each of a and b.
  600.     Clearly if ab = -q/3 is true then (wa)(wwb) == (wb)(wwa) == -q/3
  601. */
  602.  
  603.   a:=(-r/2 + Sqrt(q^3/27+ r*r/4))^(1/3);
  604.   b:=(-r/2 - Sqrt(q^3/27+ r*r/4))^(1/3);
  605.  
  606.   {a+b+adjust,w*a+ww*b+adjust,ww*a+w*b+adjust};
  607. ];
  608.  
  609. Rule("PSolve",1,1,IsUniVar(uni) And Degree(uni) = 4 )
  610. [
  611.     Local(coef4,a1,a2,a3,a4,y,y1,z);
  612.  
  613.     coef4:=Coef(uni,4);
  614.     a1:=Coef(uni,3)/coef4;
  615.     a2:=Coef(uni,2)/coef4;
  616.     a3:=Coef(uni,1)/coef4;
  617.     a4:=Coef(uni,0)/coef4;
  618.  
  619.     y1:=Head(PSolve(y^3-a2*y^2+(a1*a3-4*a4)*y+(4*a2*a4-a3^2-a1^2*a4),y));
  620.     
  621.     Concat(PSolve(z^2+1/2*(a1+Sqrt(a1^2-4*a2+4*y1))*z+1/2*(y1-Sqrt(y1^2-4*a4)),z),
  622.      PSolve(z^2+1/2*(a1-Sqrt(a1^2-4*a2+4*y1))*z+1/2*(y1+Sqrt(y1^2-4*a4)),z));
  623. ];
  624.  
  625. Function("PSolve",{uni,var})
  626.     [
  627.      PSolve(MakeUni(uni,var));
  628.      ];
  629.  
  630.  
  631. /* Generate a random polynomial */
  632.  
  633.  
  634. RandomPoly(_var,_degree,_coefmin,_coefmax) <--
  635.   NormalForm(UniVariate(var,0,RandomIntegerVector(degree+1,coefmin,coefmax)));
  636.  
  637.  
  638. /* CanBeUni returns whether the function can be converted to a
  639.  * univariate, with respect to a variable.
  640.  */
  641. Function("CanBeUni",{expression}) CanBeUni(expression,VarList(expression));
  642.  
  643.  
  644. /* Accepting an expression as being convertable to univariate */
  645.  
  646. /* Dealing wiht a list of variables. The poly should be expandable
  647.  * to each of these variables (smells like tail recursion)
  648.  */
  649. 10 # CanBeUni(_expression,{}) <-- True;
  650. 20 # CanBeUni(_expression,var_IsList) <--
  651.     CanBeUni(expression,Head(var)) And CanBeUni(expression,Tail(var));
  652.  
  653. /* Atom can always be a polynom to any variable */
  654. 30 # CanBeUni(expression_IsAtom,_var) <-- True;
  655. 35 # CanBeUni(_expression,_var)_IsFreeOf(expression,var) <--
  656.     True;
  657.  
  658. /* Other patterns supported. */
  659. 40 # CanBeUni(_x + _y,_var) <-- CanBeUni(x,var) And CanBeUni(y,var);
  660. 40 # CanBeUni(_x - _y,_var) <-- CanBeUni(x,var) And CanBeUni(y,var);
  661. 40 # CanBeUni(   + _y,_var) <-- CanBeUni(y,var);
  662. 40 # CanBeUni(   - _y,_var) <-- CanBeUni(y,var);
  663. 40 # CanBeUni(_x * _y,_var) <-- CanBeUni(x,var) And CanBeUni(y,var);
  664. 40 # CanBeUni(_x / _y,_var) <-- CanBeUni(x,var) And IsFreeOf(y,var);
  665. /* Special case again: raising powers */
  666. 40 # CanBeUni(_x ^ y_IsInteger,_var)_CanBeUni(x,var) <-- True;
  667. 41 # CanBeUni(_x ^ _y,_var)_(IsFreeOf(x,var) And IsFreeOf(y,var)) <-- True;
  668. 50 # CanBeUni(UniVariate(_var,_first,_coefs),_var) <-- True;
  669. 1000 # CanBeUni(_f,_var)_(Not(IsFreeOf(f,var))) <-- False;
  670. 1001 # CanBeUni(_f,_var) <-- True;
  671.  
  672.  
  673.  
  674.  
  675. 10 # Content(UniVariate(_var,_first,_coefs)) <-- Gcd(coefs)*var^first;
  676. 20 # Content(poly_CanBeUni) <-- NormalForm(Content(MakeUni(poly)));
  677.  
  678. 10 # PrimitivePart(UniVariate(_var,_first,_coefs)) <--
  679.     UniVariate(var,0,coefs/Gcd(coefs));
  680. 20 # PrimitivePart(poly_CanBeUni) <-- NormalForm(PrimitivePart(MakeUni(poly)));
  681.  
  682. 10 # Monic(UniVariate(_var,_first,_coefs)) <--
  683. [
  684.   DropEndZeroes(coefs);
  685.   UniVariate(var,first,coefs/coefs[Length(coefs)]);
  686. ];
  687. 20 # Monic(poly_CanBeUni) <-- NormalForm(Monic(MakeUni(poly)));
  688.  
  689. 30 # Monic(_poly,_var)_CanBeUni(poly,var) <-- NormalForm(Monic(MakeUni(poly,var)));
  690.  
  691.  
  692. 10 # BigOh(UniVariate(_var,_first,_coefs),_var,_degree) <--
  693.     [
  694.      While(first+Length(coefs)>=(degree+1) And Length(coefs)>0) DestructiveDelete(coefs,Length(coefs));
  695.      UniVariate(var,first,coefs);
  696.     ];
  697. 20 # BigOh(_uv,_var,_degree)_CanBeUni(uv,var) <-- NormalForm(BigOh(MakeUni(uv,var),var,degree));
  698.  
  699.  
  700.  
  701. Horner(_e,_v) <--
  702. [
  703.   Local(uni,coefs,result);
  704.   uni := MakeUni(e,v);
  705.   coefs:=DestructiveReverse(uni[3]);
  706.   result:=0;
  707.  
  708.   While(coefs != {})
  709.   [
  710.     result := result*v;
  711.     result := result+Head(coefs);
  712.     coefs  := Tail(coefs);
  713.   ];
  714.   result:=result*v^uni[2];
  715.   result;
  716. ];
  717.  
  718.  
  719. DivPoly(_A,_B,_var,_deg) <--
  720. [
  721.   Local(a,b,c,i,j,denom);
  722.   b:=MakeUni(B,var);
  723.   denom:=Coef(b,0);
  724.  
  725.   if (denom = 0)
  726.   [
  727.     Local(f);
  728.     f:=Content(b);
  729.     b:=PrimitivePart(b);
  730.     A:=Simplify(A/f);
  731.     denom:=Coef(b,0);
  732.   ];
  733.   a:=MakeUni(A,var);
  734.  
  735.   c:=FillList(0,deg+1);
  736.   For(i:=0,i<=deg,i++)
  737.   [
  738.     Local(sum,j);
  739.     sum:=0;
  740.     For(j:=0,j<i,j++)
  741.     [
  742.       sum := sum + c[j+1]*Coef(b,i-j);
  743.     ];
  744.     c[i+1] := (Coef(a,i)-sum) / denom;
  745.   ];
  746.   NormalForm(UniVariate(var,0,c));
  747. ];
  748.  
  749.